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Abstract 



This paper is concerned with the numerical solution of nonlinear ill-posed op- 
erator equations involving convex constraints. We study a Newton-type method 
which consists in applying linear Tikhonov regularization with convex constraints 
to the Newton equations in each iteration step. Convergence of this iterative reg- 
ularization method is analyzed if both the operator and the right hand side are 
given with errors and all error levels tend to zero. Our study has been mo- 
tivated by the joint estimation of object and phase in 4Pi microscopy, which 
leads to a semi-blind deconvolution problem with nonnegativity constraints. The 
performance of the proposed algorithm is illustrated both for simulated and for 
three-dimensional experimental data. 

1 Introduction 

In this paper we present and analyze a Newton-type regularization method for 
nonlinear ill-posed operator equations with convex constraints. More specifically, 
let X and y be Hilbert spaces, C C X a nonempty, closed convex set, and 
F : C ^ y a "forward" operator, which we assume to be Gateaux differentiable. 
We consider the inverse problem of reconstructing in the operator equation 



if only noisy versions of both F and g are given. Moreover, we aim to prove 
convergence of such reconstructions as the noise levels tend to zero. 

An inverse problem for which it is particularly important to properly incor- 
porate a convex constraint into the inversion scheme arises in a confocal fluores- 
cence microscopy technique (cf. |17j ) called 4Pi microscopy. This technique was 
suggested and developed by Hell et.al. [91 [10] and allows for a substantial en- 
hancement of resolution using interference of two laser beams in the microscopic 
focus and/or interference of fluorescence photons on the detector. In standard 
confocal microscopy the relation between the unknown fluorescent marker density 




x 




(1.1) 



1 



/ € L 2 (M 3 ) of the specimen and the measured intensity g is given by a convolu- 
tion with a point spread function (psf) h G L 1 (M 3 ), which is often modeled as a 
Gaussian function: 

S(x)= /" fc(x- y )/(y)dy (1.2) 

The width of /i is typically much larger along the so-called optical axis (which we 
assume to be the X3-axis) than in directions perpendicular to the optical axis. 

4Pi microscopy allows an increase of resolution along the optical axis by a fac- 
tor of 3-7 using interference of coherent photons through two opposing objective 
lenses. Here the psf is no longer spatially invariant in general, but depends on the 
relative phase </>(x) of the interfering photons, which has to be recovered together 
with the fluorophore density / in general since it depends on the refractive index 
of the specimen which is unknown. The imaging process can be modeled by an 
operator equation i*4Pi(/, (ft) = g with a forward operator of the form 

F 4Pi (/, 0)(x) := J p(y - x, <^(x))/(y)dy . (1.3) 

Note that F is nonlinear in (ft and that / i— > F (f,4>) is not a convolution operator 
in general. As a density, / has to be nonnegative. Therefore, we have the convex 
constraint (/, (ft) G C with C := {(/, (ft) : / > 0}. A simple frequently used model 
for the 4Pi-psf (cf., e.g., [I]) is given by 

p(x, if) « fc(x) cos n (cx 3 + I) , (1.4) 

where h is the psf of the corresponding confocal microscope, and the cosine term 
represents the interference pattern for different types of 4Pi-microscopes corre- 



sponding to n = 2,4, respectively (see Fig. 1.1). So far reconstruction of / in 
commercially available 4Pi microscopes is done by standard deconvolution soft- 
ware assuming the relative phase function (ft to be constant. Although spatial 
variations of (ft can approximately be avoided experimentally in some situations, 
the assumption that (ft is constant imposes severe limitations on the applicability 
and reliability of 4Pi microscopy. Therefore, it is of great interest to develop 
algorithms for the solution of the convexly constrained nonlinear inverse problem 
to recover both the object function / and the relative phase function (ft from the 
data g. 

To this end we propose and analyze the following constrained version of the 
iteratively regularized Gaufi-Newton method (IRGNM). We assume that both 
the right hand side g in the operator equation ( |1.1[ ) and the operator F are only 
given approximately with errors by g$ and F$, respectively. Error bounds will be 
specified in the next section. Given some initial guess xo G C, we consider the 
iteration 



%n+i = argmm 



\F$[x n ](x - x n ) + F s (x n ) - gs\\ 2 + a n \\x - x \\ 2 , (1.5a) 



2 



••• 



•I* •••• 



(a) if = 



II* 



(d) cp = 



(b)p 



(e) 99 



(c) 93 



TT 



7T 



Figure 1.1: The top line shows the psf of a 4Pi microscope modeled by (1.4) for relative phases 
(p = 0, 5,7r on the plane containing the optical axis, which is indicated by the white arrow. 



The bottom line shows the more accurate model (3.1). 



n = 0, 1, 2, . . . , with a sequence of regularization parameters a n satisfying 

Oi 

1 < < r, lim a n = 0, a n > for some r > 1 and for all n £ No- 

a n+ i n-+oc 

(1.5b) 

In the unconstrained case C = X this reduces to the IRGNM as suggested in the 
original paper by Bakushinskii [2j . For C 7^ X a quadratic minimization problem 
with convex constraint has to be solved in each Newton step. In [2J convergence 
rates were shown for Holder type source conditions with exponent v = 1. In [Hll2j 
order optimal convergence rates for more general Holder type and logarithmic 
source conditions were proven. For numerous further references on the IRGNM 
and other iterative regularization methods we refer to the monographs [31 115] . 
More recently, Kaltenbacher and Hofmann [H] proved optimal convergence rates 
of the IRGNM in Banach spaces for general source conditions. 



The convergence result we will present in the next section (Theorem 2.1) 
takes into account two features, which are essential for 4Pi reconstructions and 
are not covered in the literature so far: First of all, our source condition takes 
into account the convex constraint and is weaker than the corresponding source 
condition for the unconstrained case, yielding the same rate of convergence. This 
reflects the observation reported below that projecting reconstructions of the 
unconstrained IRGNM onto C does not yield competitive results. For linear 
Tikhonov regularization with convex constraints we refer to Neubauer [16] and [0, 
section 5.4]. Moreover, unlike many other references on the IRGNM, we also take 
into account errors in the operator since they are important in our application: 



The frequently used model (1.4) for the 4Pi psf is only a first approximation, and 
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even the more accurate model based on the evaluation of diffraction integrals, 
which we used in our code (see Fig. |l.l| and eq. (3.1 ) below), contains parameters, 
which have to be estimated including errors. Other references discussing the 
influence of errors in the operator for the IRGNM include [3J and |13j . 

The plan of this paper is as follows: Our main convergence result, Theorem 
2.1 is formulated and proved in Section[2j Section [3] contains a more detailed dis- 



cussion of 4Pi microscopy and the model ( |1.3[ ), a comparison with other methods, 
and numerical results both for simulated and experimental data. 



2 IRGNM with Convex Constraints 
2.1 Formulation of the theorem 

We assume that F, F$ : C — > y are both Gateaux differentiable with bounded 
derivatives Fs[x] for all x G C and that the following error bounds hold: 



\\g - gs\\ < 5 g , 
F(x j )-F 5 (x j ) <5 F 

F'[x^]-F' 5 [x^] <5 F , 



with noise levels <5 5 ,<5f,5f/ > 0. 

Further we assume that a source condition of the form 
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Pc{F' [x'\* uj + xq) for some to £ y with ||w|| < p 



(2.1a) 
(2.1b) 

(2.1c) 
(2.2a) 



is satisfied where Pq '■ X — > C denotes the metric projection onto C. The source 
condition (2.2a) corresponds to t he one for linear constrained Tikhonov regular- 

n((T*T)V 2 ) 



2.2 



and since 1Z(T*) 



for a bounded 



ization we assume in Lemma 
linear operator T : X — » y (cf. |6., Proposition 2.18]) it corresponds to a Holder- 
type source condition with exponent v = \. As the ( |2.2a ) contains the projector 
Pc, it is less restrictive than in the unconstrained case C = X . In particular, x^ 
may not be smooth even if i^'fa^]* is smoothing and xq is smooth. 
If -F'[x'] is not injective, we further assume that x^ satisfies 



x' = argmm ||x 

{a;6C:F'[a;t](a;-2:t)=0} 



(2.2b) 



Obviously, this condition is empty if F'^] is injective. Moreover, if for any 
vq G N{F'[x^\) there exists a differentiable curve v : [0,e) — > C with v(0) = 



x\ v'(0) = vq and F(v(t)) = g for all t (see e.g. [8] for a problem where 
this condition is satisfied), then it is easy to see that (2.2b) follows from x^ = 
argmin {:r6C:F{:E)=9} ||x - x \\. 

As nonlinearity condition on the operator Fg we only need to assume that for 
some 7 > there exists a Lipschitz constant L > such that 



-tl 



< L 


x — x^ 







for all x G C with 



< 7- (2-3) 



We can now formulate our main convergence theorem: 
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Theorem 2.1. Assume that (1.1) and ( 2.1 ) — ( 2.3 ) are satisfied with p is suffi- 
ciently small, set 



5 := max(5 g + <$p, <5fv 



(2.4) 



and consider the sequence (x n ) defined by (1.5). 

Then the iterates satisfy \\x n — ar < 7 and in the noise free case 5 = we have 



Xn X 



n — > 00. 



For 5 > assume that a stopping index N is chosen such that 

ct 7v < r]5 < a n , < n < N 



(2.5) 



(2.6) 



with some constant r/ > sufficiently large. Then the error of the final approxi- 
mation fulfills 



xn — x 



o vs 



5^0. 



(2.7) 



2.2 Proof of the theorem 



Note that if F = F$ = T : X — > y is linear and bounded, then (1.5a) reduces to 
linear constrained Tikhonov regularization 



x a := argmm 

x&C 



\Tx — gs\\ 2 + a \\x — xq\\ 2 



for a sequence of regularization parameter a = a n . We first recall the stability 
and approximation properties in this case since they will be needed later in the 
proof. 



Lemma 2.2. 1. If 



x a := argmm 



\Tx — g s \\ 2 + a \\x — xq\\ 2 



for some g s G y , then 



I *^ Qt *^ OL 1 1 



\98 ~ 9s\ 
a 



(2.9) 



2. Let g = g s G T{C) and x$ G C, and assume that the best- approximate- 



solution x c := argmin| 2 . eC . r;r=9 } [|x — x§\ satisfies the source condition 



for some uj G y . Then 



Xrv X( 



4 = Pc(T*uj + s ) 



\/a||w|| and \\Tx a — g\\ = a \\uj\ 



(2.10) 



(2.11) 
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Proof. In (2.9) the special case xo = is proved in [6j Theorem 5.16]. The general 
case can be reduced to this special case by the substitution of variables z = x — xq 
since 



Xa — %o = argmm 

z&C—xq 



\Tz + Txq — gs\\ 2 + a \\z\\ 2 



(2.12) 



(2.11) can be reduced to the case xq = 0, which is covered by [6] Theorem 5.19], 
by the same substitution of variables and the identity Pc_ X0 (T*uj) = Pc(T*oj + 

Xo) - Xq. □ 

Next we need a stability estimate with respect to perturbations of the opera- 
tors, i.e. an estimate on the difference of 



argmm 



\Ti(x — x)\\ 2 + a \\x — xq\\ 2 



, i€{l,2}, (2.13) 

where T\,T2 : X — > y are bounded linear operators and a > 0. Using the 
optimality conditions for the minimizers of (2.13), a straightforward computation 
gives an estimate of the form 



\xi - x 2 \\ < - \\Ti - T 2 \\ . 
a 



(2.14) 



This simple estimate is not sufficient for our purposes, however. The follow- 
ing proposition shows that under a source condition we can obtain an improved 
estimate with a constant independent of a: 



Proposition 2.3. Let x\ and x 2 be defined by (2.13). Moreover, let the source 
condition 



x = PciT^uj + x ) 



(2.15) 



hold for some uj € y and let x = argmin| :rgC:T2;c=r2j .| — xq\\. Then the distance 
of xi and x 2 is bounded by 



\xi - x 2 \\ < 



T 2 \ 



Proof. Prom Lemma 2.2, part 2 we obtain 

||T 2 (sc2 — x)\\ < a\\u\\ 
\\x 2 — x\\ < \fa ||u;|| . 

Let xc '■ X — > M + be the proper, convex and lower semicontinuous functional 



(2.16a) 
(2.16b) 



Xc{x) :-- 



x e C 
oo otherwise 



and let qi £ dxc(xi), then for i E {1,2} the first order optimality condition for 
the minimizers Xi is given by 



T*Ti(xi - x) + a(xi - xo) + qi = 0. 



(2.17) 
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Subtracting the equations (2.17) gives 



r*Ti(xi - x 2 ) + a(xi - x 2 ) + (qi - q 2 ) = (T 2 T 2 - T{Ti)(x 2 - x). 

Now taking the inner product with x\ — x 2 we obtain 

\\T\{xi - x 2 )\\ 2 + a \\xi - x 2 \\ 2 + (qi - q 2 ,x 1 - x 2 ) = 

((T 2 * - Tl)T 2 (x 2 - x),x\ - x 2 ) + ((T 2 - Ti)(x 2 - x),T 1 (xi - x 2 )). 

The right hand side can be estimated with help of Young's inequality 

((T 2 * - T*)T 2 (x 2 - x),xi - x 2 ) + ((T 2 - T 1 )(x 2 - x),Tx{x\ - x 2 )) 



(2.18) 



(2.19) 



< — ||Ti - T 2 
~ 2a 11 



| 2 \\T 2 (x 2 - x)\\ 2 + ^ \\xi - x 2 \\ 2 + - ||Ti - T 2 \\ 2 \\x 2 - x\\ 2 + \\Ti(xi 



Using (qi — q 2 ,x\ — x 2 ) > (see, e.g., Section 9.6.1, Theorem 1]), (2.19) and 
the inequalities (2.16) the assertion follows from 

77 || - x 2 || 2 < — \\T 2 (x 2 - x)\\ 2 \\Ti - T 2 \\ 2 + - \x 2 - 5|| 2 ||Ti - T 2 || 2 
I la 1 

< ^a \\u\\ 2 \\Ti - T 2 \\ 2 . 

□ 

Now we are able to formulate a recursive error estimate for the IRGNM with 
closed convex constraint. 



Lemma 2.4. Assume that (1.1) and (1.5)-(2.3) are satisfied and that x n G C 
with \\x n — x^\\ < 7, then the error e n := x n — x^ satisfies 

1 L Hi 1 [3 

1 1 6n+i || — 1 — 77 1 1 

|| + \ -pL \\e n \\ + —— (5 g + 5 F ) + y/OnP + \ -p§F'- 



foun 2 



'a r 



(2.20) 



Proof. At first we note that one can express the noisy data as gs = Fs(x^) +£ + e, 
with ||£ || < 5f and ||e|| < 5 g . Further since x n S C and F$ is Gateaux differentiable 
with derivatives that fulfill condition ( 2.3 ), we can express Fs(x^) in a Taylor series 

F 5 (xt) = F & {x n ) + F' s [x n ]{x^ -x n )+r(x^ -x n ), (2.21) 

where 

L 



r(x^ — x n ) 



< 



(2.22) 



Thus we can rewrite the IRGNM functional (1.5a) of the ra-th iteration step, 
defining T n := F' s [x n ], as 



\T n x - (T n x n - F s (x n ) + gs)\\ 2 + a n \\x - xq\ 



T n (x - 2+) - r(x t - x n ) - £ 



+ a n \\x — xq\ 



(2.23) 
(2.24) 



7 



Now we can decompose the distance of the solution x n +i of the (n+l)-th iteration 
to the exact solution x* using the triangle inequality 



with 



x n+ 1 = argmm 

x£C 



argmm 

x&C 



T n (x - x*) - r n (x^ - x n ) - £ - e 



+ a n \\x — x \ 



(2.25) 



T n (x - x ] ) 



+ a n \\x — xq\ 



argmm 

x&C 



F'[x^]{x-x^) 



+ a n \\x — xq\ 



It follows from Lemma 2.2, part 1 that 



\r n (x^ - x n ) + £ + el 



'OLr. 



With (2.22), ||f || < 5 F and ||e|| < 5 g we obtain 



l^n+l ^o n ,n|| 5: 



1 ( L 



lou,. V 2 



X n X 



(2.27) 



(2.28) 



The second term in (2.25) can be estimated using Proposition 2.3 with T\ = T n , 
T2 = F'[x^] and x = x\ which gives 



T n -F' s [x^ + F' 5 [x^-F'[x^ 
+ S F 



(2.29) 



where we used (2.3) and (2.1c) to obtain the last inequality of (2.29). For the 
third term in (2.25) we again use Lemma 2.2, part 2 to obtain 



X fy„ X 



< \J~Un~P- 



Combining (2.25), (2.28), (2.29) and (2.30) gives the assertion. 



(2.30) 

□ 



Estimate (2.20) is of the form used in [3] and [TJ], so we can use a similar 
proof now to obtain the main result. 



Proof of Theorem 2.1, It follows from Lemma 
fulfill the inequality 



2.4 



that the quantities 0„ := 



9 n+ i < a + bG n + cQ 2 n 



(2.31) 
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with a := y/rp for 5 = and a := \/r(p+ 



Sg+Sp 



+ 



y-^fL) for 5 > 0, b 



%rpL 



and c := s/r-^- Let ti and £2 be solutions to the fixed point equation a+bt+ct = t, 
i.e. 



2a 



1-6+7(1^ 



4ac 



t2= 1 - tW( l- tF -4 ^ (232) 



let the stopping index iV < 00 be given by (2.6) and define Cq := max(0o,ii). 
We will show by induction that 



e n <c© 



for < n < N if 



b + 2y/ac < 1, 
G < i 2 , 
\/aoCe < 7. 



(2.33) 



(2.34a) 
(2.34b) 
(2.34c) 



Conditions (2.34) are satisfied if p is sufficiently small and r\ sufficiently large. 



For n = (2.33) is true by the definition of Cq. Assume that (2.33) is true for 



some k < N, then by (2.34c) and Lemma 2.4, (2.31) is true for n = k. Condition 



(2.34a) assures that ii,i2 £ K and t\ < t2, and by (2.33) one has < 0^ < t\ or 



ti < 0fc < ©o- In the first case, since a,b,c> 0, we obtain 

©fe+i < a + bG k + cG^ < a + fai + ct? = t x . (2.35) 

l)i + ct 2 < for ti < t < t 2 



In the second case by (2.34b) and the fact that a + 
we obtain 



6 fc+ i < a + b& k + c&l <9 k < G . 



(2.36) 



Thus in both cases (2.33) holds for n = k + 1 and the induction is complete. 



By definition N = 00 for 6 = and thus (|2.33|) directly implies (2.5). Using 

□ 



O-N < by (2.6) also assertion (2.7) follows directly from (2.33) 



3 Joint reconstruction of object and phase 
in 4Pi-microscopy 

Before we describe our mathematical model of the 4Pi imaging process precisely, 
let us discuss this technique in some more detail. Confocal fluorescence mi- 
croscopy allows the reconstruction of three-dimensional fluorescent marker den- 
sities in living cells by scanning a specimen at a grid of points {xj £ I 3 : j = 
1, . . . , iV}. Laser light is focused to a small area by objective lenses, and a pinhole 
is used to collect only fluorescence photons emitted close to the focus xj (cf. p2] ) . 
The psf /i(x — y) in (1.2) is the probability that a fluorescence photon emitted 
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at y is detected if the point x is illuminated. Data consist of photon count num- 
bers Gj, j = 1, . . . , N, which are Poisson distributed random numbers with mean 

In 4Pi microscopy the same data model holds true, but g is given by g = 



i*4Pi(/, 4>) with the integral operator (1.3). For its kernel we use the more accurate 
model 

p(z, <p) = \Ex(z) - exp(i0)E 2 (z)| 2 /i det (z) (3.1) 



(see [IB]) instead of the simple model (1.4). Here Ej^ are counterpropagating 



focal fields and fi^et is the detection psf, for which we used implementations 
available under www . imspector . de. 

3.1 forward operator and its derivative 



We first define appropriate function spaces for the integral operator F$p\ in (1.3). 
We assume that / is supported in some cube Vt := Y\^ = i[—Rj, Rj] and choose 
L 2 (Q) with the standard L 2 -norm as function space for the object /. We may 
further assume that p(-,f) is supported in some (typically much smaller) cube 
n^it - r ji r j] f° r an V 9 such that g is supported in £1' := nj/=i[ — Rj ~ r j'Rj + r j\- 
A reason why joint reconstruction of / and (f) from data g often works even 
though the problem is formally underdetermined, is that 4> can be assumed to be 
very smooth (often it is even assumed to be constant). Therefore we choose the 

Sobolev space H 2 {VL') for <p with norm \\4>\\ H 2 {nl) := (f n , \(p\ 2 + |V>| 2 + \ A(p\ 2 dx)* 
to achieve smooth interpolation in areas where no information on <fi is contained 
in the data. (This is the case, e.g., in areas where / is constant. But in such 
areas 4> is irrelevant for the primary goal to recover /.) 

The data misfit term should reflect the distribution of data errors. Since our 
data are Poisson distributed, a natural data misfit term would be the negative 
log-likelihood function, which is given by YljLi 9( x j) — Gj l g5( x i)- We define a 
piecewise constant approximation g s € L 2 {£1') of the data (Gj) and approximate 
the negative log- likelihood by a second order Taylor expansion at Gj. This leads 
to a weighted L 2 space y := L 2 (Q',w) with norm ||<7||y = f n , g 2 (x)w(x) dx and 
weight function 

v ' 2max(# 5 (x),c) 
Here c > is a small constant avoiding division by zero. As usually multiple 
weaker sources contribute to the data noise, a suitable choice of c is the back- 
ground noise level. Better approximations to the Poisson log-likelihood can be 
achieved by taking a Taylor expansion at g n = i ? 4Pi(/ n) ^n) and iterating in a 
sequential quadratic programming manner, but for the count rates in our exper- 
imental data this did not lead to a noticible improvement. 

In summary, the precise definition of our forward operator is as follows: 

F 4Pi : L 2 (Q) x H 2 (n') — ► L 2 (Q',w) 

(F 4Pi (/,0))(x) := /p(x-y,0(x))/(y)dy, x € O' . ( "' 2) 
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Note that i*4p; does not change if p(-, ip) is replaced by its periodic extension with 
period cell £2' . 

Lemma 3.1. If p : Y\ 3 j=i( R /( 2 ( R j + x (K/ttZ) -)■ M is continuous and 

continuously differentiable with respect to its last argument, then the operator 
F^pi defined in (3.2) is Frechet differentiable on X with 

*4Pi[/,0(fc/,fy)( x ) = i £{Ky-x,0(x))^ / (y) + ^(y-x,0(x))/(y)^(x)jdy, 

(3.3) 

and £/ie adjoint of F'[f,(f>] : L 2 (f2) x H 2 {Q.) — )> L 2 (fy,u;) is giwen &?/ 

/ U(-xi(x))j(x) W (x)«ix \ 

FW/ ' 0] 9= { f ( 9W f n d 4(--yM-))f(y)dy) ) (3 ' 4) 

where j : H 2 (Q') L 2 ($7') is i/ie embedding operator. Moreover, F^ Fi satisfies 
the Lipschitz condition (2.3) i/ ^| is uniformly Lipschitz continuous with respect 
to its last argument. 



Proof (sketch). The Frechet differentiability of F and eq. (3.3) follow from a 
Taylor expansion of the kernel p with standard estimates on the Taylor remainder 
and the continuity of the embedding H 2 (9!) L°°(9!). The adjoint F'[f, ^]* 2 of 
the continuous extension F'[f, 4>] L 2 of F'[f, eft] to L 2 (f2) x L 2 (S7') can be computed 
by interchanging the order of integration. Then (3.4) follows from F'[f,<f)]* = 
(F'[f,4>] L 2 (o°j))* = (o j*) F'ifi^LZ- Tri e statement on Lipschitz continuity is 
straightforward. □ 

The crucial observation for an efficient implementation of F$p\ and F^ Pi is 
that p can be separated into 

M 

p(z,cp)= ex?p{im<p)A m (z) 

m=-M 

with A m G L 2 (n|=i(K/(2( J R J ' + 7j)Z))). This was observed by Baddeley et al. 
in pQ for the approximation (1.4) and by Vicidomini et al. in [18] for the model 
(3.1). Hence, 

M 

(i ? 4Pi(/,0))(x)= ]T exp(im^»(x)) / A w (x-y)/(y)dy, xefl'. 

m=-M Jn ' 

Here / is extended by in fi' \ (zero-padding). The convolution integrals can 
be evaluated efficiently using FFT. An analogous procedure can be applied for 
the evaluation of F'[f,(p] and its adjoint. 

We approximated the phase <j> using tensor products of Chebychev polynomi- 
als, for which the Gramian matrix with respect to the H 2 inner product can be 
computed explicitly. 
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3.2 Implementation and necessity of the nonnegativ- 
ity constraint 

We solve the constrained quadratic minimization problems 

(fn+i,4>n+i) ■= argmin || F' iVi [(f n , (f> n )](f, 4>) - g n \^ + a n \\(f, <p) - (fo,<po)\\ 2 , 
(f,it>)&L 2 (n)xH 2 (n), 

/>0 a.e. 

(3.5) 



with 



9n ■= K Pi [(f n ,(p n )}((f n ,(j) n )) - F 4Pi ((f n ,(p n )) + g s 



using the semi-smooth Newton method (cf. [11]). In each step of this method 
an unconstrained, positive definite linear system has to be solved, which is done 



by the conjugate gradient method. In Figure 3.1 the reconstruction of a fluo- 



rophore density and the phase from a 2d-slice of real 4Pi data is shown. To 




(c) object reconstruction /12 (d) phase reconstruction $12 



Figure 3.1: Panel (a) shows a slice of real 3-dimensional 4Pi-data, where in the approximate 
center the waves were interfering destructively. From these data the reconstructions of object 
|(c) I and phase [(d)] have been obtained with the constrained IRGNM for /o = and </>o = 0. 
The reconstruction of the phase reflects the constructive interference on the left and right side 
of the data and the destructive interference in the center. The modeled psf (for constructive 
interference) is depicted in panel |(b)| together with indications on the scale and the optical 
axis (represented by the arrow). To reconstruct the phase we used a basis of polynomials 
with maximal degree 7 in each dimension. 
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Figure 3.2: Panels (a) and (b) show reconstructions of object and phase respectively from the 
data shown in Figure |3.1a[ which have been obtained by the unconstrained IRGNM. The 
phase is very badly reconstructed which leads to remaining sidelobes in the object reconstruc- 
tion. This is most eminent in the center where the psf features destructive interference. The 
same deficiencies can be observed in panels (c) and |(d)[ where the nonnegativity constraint 
has been incorporated by a simple projection after each step. All of the reconstructions 
depicted in Figures 3.1 and 3.2 were performed with the same regularization parameters. 



demonstrate the necessity to incorporate nonnegativity constraint into the mini- 



mization problem, in Figure 3.2 we display reconstructions from the same data, 
without constraint and with simply projecting onto C after an unconstrained 
IRGNM step, i.e. for the iteration schemes 

1 1 2 

(/n+nC+i) : = argmin M|i^ 4 Pi[(/^ ,<0(/, 4>) ~ 9n\\ +"«!!(/, 0) - (fo,<Po\ 
(fn+v€+i) ■= ^cargmin [||Fl Pi [(/P, <0(/, 0) " 9n\\ 2 + a n \\(f,<f>) - (/ o ,0o)|| 2 
Here the metric projection is given by Pc(f,4>) = (max(/, 0), <ft). Comparing the 



reconstructions of Figures 3.1 and |3.2| it is obvious that the incorporation of the 
nonnegativity constraint in the minimization problem is necessary for accurate 
reconstructions of the phase. 

A further option pursued in [18J is to update / and <j> in alternating manner 
such that in each update step for / a constrained minimization problem for / only 
instead of both / and <f> has to be solved. However, such a procedure requires 
significantly more iteration steps. 
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3.3 Results for simulated and experimental data 




(a) simulated object /' (b) simulated phase 



(c) noisy data g$ 



(d) object reconstruc- 
tion /27 for noisy data 




(e) phase reconstruc- 
tion 027 for noisy data 



* residual 

* object error 

* phase error 
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(f) errors and residual for 
noisy data 




(g) object reconstruc- 
tion for exact data 



(h) phase reconstruc- 
tion </>3fi for exact data 



* residual 

* object error 

* phase error 



5 10 15 20 25 30 35 

(i) errors and residual for ex- 
act data 



Figure 3.3: Panels (a) and (b) show a simulated object and phase, and panel (c) the corre- 
sponding 4Pi-data perturbed with Poisson-noise. Panels [(d)] and [(e)] show the reconstructions 
of object and phase respectively from this noisy data obtained with the constrained IRGNM. 
In panel (f) the residual \\gs — F(f n ,(/) n )\\ L 2, the object error |L/n~ J_ an d the phase error 



L 2 are plotted over the iteration index n. Panels (g) - (i) are analogous to panels 



(d) - (f) with noisy data g$ replaced by exact data F(f', <jr). 



Figure 3.3 shows reconstructions from simulated two-dimensional noisy and 
exact data. Here we chose polynomials of maximal degree 7 in each dimension to 
approximate the reconstructed phase. We chose the exact phase as a shifted sum 
of a sine and arctan function, which does not belong to the polynomial subspace. 
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For exact data the sidelobes are removed completely, and for noisy data at the 
given count rate only very little of the sidelobes is left in the reconstruction. The 
required number of semi-smooth Newton (SSN) steps increases with n. To give 
an idea, we mention that less than 8 SSN steps were needed for n < 21 with less 
than 80 CG steps in each SSN step, and for n = 30 the algorithm required 49 
SSN steps with less than 600 CG steps. We must say that the stopping indices 
for the Gaufi-Newton iteration are chosen somewhat arbitrarily in this paper. 
The development of a good stopping rule for the kind of errors considered in 
this paper, nonlinear operators and convex constraints is an interesting topic for 
future research. 




(a) simulated object /' (b) simulated phase 



■ « » i $t : i'ii - 



(c) noisy data g$ 
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(d) object reconstruc- 
tion /i4 



(e) phase reconstruc- 
tion 014 



(f) EM-TV object re- 
construction 



Figure 3.4: Panels (a) and (b) show a simulated object and phase from which noisy 4Pi 
data was created (panel (c)[). Panels |(d)| depicts the object reconstruction obtained with 
the IRGNM and panel (e) the corresponding phase reconstruction. Panel [(f)] shows a recon- 
struction by an expectation maximization algorithm with TV penalty using the reconstructed 
phase (i 



In Figure 3.4 we chose an object which is constant in a region, and hence the 
data carry no information on the phase there. Due to the -£f 2 -phase penalty term, 
the phase is interpolated smoothly in this area and recovered quite well, except 
in dark areas close to the boundary. In contrast, the reconstruction of the object 
exhibits a grainy structure in the central area. This is a consequence of choosing 
the I? norm as object penalty. Since we have found a good approximation </> ap p 
of the phase, we can compute a better reconstruction of the object in a second 
step by solving an inverse problem for the linear operator / i— > F(f, app ). The 
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result in Figure 3.4e was computed using an expectation- maximization method 
with a TV penalty term and Bregman iterations as described in [5]. 

Figure |3,5| shows cuts through 3-dimensional experimental data. The corre- 



sponding reconstructions of object and phase are shown in Figure 3.6. Note that 
due to the simultaneous reconstruction of the phase function the non-symmetric 
sidelobes in the data have been removed in the object reconstruction. Including 





Figure 3.5: Data of microtubules in a Vero cell, for NA = 1.34, A cx = 635nm, A cm = 680nm. 
The data extension is (2952nm x 9296nm x 1904nm) in x, y and z direction respectively. The 
annotations at the axes number the voxels in the respective dimension. 

the zero padding, the data contained approximately 2 million voxels. The phase 
has been approximated by Chebychev polynomials of maximal degree 3 in each 
dimension. 
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(a) object reconstruction 



(b) phase reconstruction 



Figure 3.6: Panel 



3.5 In Panel (b) 



(a) shows the reconstruction of the object from the data shown in Figure 
the corresponding reconstruction of the phase is depicted. 
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